1 Dərsin Məqsədləri

Bu dərsin sonunda tələbələr:

GRM mathematical formulation və polytomous IRT prinsiplərini anlaya bilə
Category Response Functions (CRFs) yarada və interpretasiya edəcək
Threshold parameters və onların mənaını izah edəcək
Rating Scale ModelPartial Credit Model ilə GRM-i müqayisə edəcək
R paketlərində (ltm, mirt, eRm) polytomous IRT analysis aparacaq
Model fitmodel selection metodlarını tətbiq edəcək
Practical applications ilə real data-da GRM istifadə edəcək

2 GRM-nin Riyazi İfadəsi

2.1 Graded Response Model Mathematical Foundation

cat("=== GRM MATHEMATICAL FOUNDATION ===\n\n")
## === GRM MATHEMATICAL FOUNDATION ===
cat("Graded Response Model: Ordered categorical responses üçün IRT model\n")
## Graded Response Model: Ordered categorical responses üçün IRT model
cat("• Samejima (1969) tərəfindən yaradılmış\n")
## • Samejima (1969) tərəfindən yaradılmış
cat("• Cumulative probability approach istifadə edir\n")
## • Cumulative probability approach istifadə edir
cat("• Discrimination və threshold parameterlər\n")
## • Discrimination və threshold parameterlər
cat("• Likert scales və rating data üçün ideal\n\n")
## • Likert scales və rating data üçün ideal
# Comprehensive GRM mathematical demonstration
demonstrate_grm_mathematics <- function() {
  
  cat("=== GRM MATHEMATICAL DEMONSTRASIYASI ===\n")
  
  # Generate example GRM data
  grm_data <- generate_grm_data(
    n_items = 6,
    n_persons = 1500,
    n_categories = 5,  # 0, 1, 2, 3, 4 (5-point Likert)
    seed = 1601
  )
  
  cat("GRM Example Dataset:\n")
  cat("Items:", ncol(grm_data$responses), "\n")
  cat("Persons:", nrow(grm_data$responses), "\n")
  cat("Categories per item:", grm_data$n_categories, "(0-4)\n")
  cat("Response range:", range(grm_data$responses), "\n\n")
  
  # Display mathematical formulation
  cat("=== GRM MATHEMATICAL FORMULATION ===\n")
  
  mathematical_formulation <- c(
    "**CUMULATIVE PROBABILITY MODEL:**",
    "",
    "P*ik(θ) = P(Xij ≥ k | θi) = exp[ai(θi - bik)] / [1 + exp[ai(θi - bik)]]",
    "",
    "Where:",
    "• P*ik(θ) = Cumulative probability of responding in category k or higher",
    "• ai = Discrimination parameter for item i",  
    "• bik = Threshold parameter k for item i",
    "• θi = Ability level of person i",
    "",
    "**CATEGORY PROBABILITIES:**",
    "",
    "P(Xij = k | θi) = P*ik(θ) - P*i(k+1)(θ)",
    "",
    "Special cases:",
    "• P(Xij = 0 | θi) = 1 - P*i1(θ)",
    "• P(Xij = mi | θi) = P*imi(θ)  [highest category]",
    "",
    "**CONSTRAINTS:**",
    "",
    "• bi1 < bi2 < ... < bi(mi) [Ordered thresholds]",
    "• Σk P(Xij = k | θi) = 1 [Probabilities sum to 1]",
    "• All P(Xij = k | θi) ≥ 0 [Non-negative probabilities]"
  )
  
  for(line in mathematical_formulation) {
    cat(line, "\n")
  }
  
  # Demonstrate with specific example
  cat("\n=== NUMERICAL EXAMPLE ===\n")
  
  # Use first item from generated data
  example_item <- 1
  a_example <- grm_data$discriminations[example_item]
  b_example <- grm_data$thresholds[example_item, ]
  
  cat("Example Item", example_item, ":\n")
  cat("Discrimination (a):", round(a_example, 3), "\n")
  cat("Thresholds (b):", paste(round(b_example, 3), collapse = ", "), "\n")
  
  # Calculate probabilities at different theta levels
  theta_examples <- c(-2, -1, 0, 1, 2)
  
  probability_table <- data.frame()
  
  for(theta in theta_examples) {
    
    # Calculate cumulative probabilities
    cum_probs <- numeric(grm_data$n_categories)
    for(k in 1:(grm_data$n_categories - 1)) {
      cum_probs[k] <- exp(a_example * (theta - b_example[k])) / 
                      (1 + exp(a_example * (theta - b_example[k])))
    }
    cum_probs[grm_data$n_categories] <- 0
    
    # Calculate category probabilities
    cat_probs <- numeric(grm_data$n_categories)
    cat_probs[1] <- 1 - cum_probs[1]
    
    for(k in 2:(grm_data$n_categories - 1)) {
      cat_probs[k] <- cum_probs[k-1] - cum_probs[k]
    }
    cat_probs[grm_data$n_categories] <- cum_probs[grm_data$n_categories - 1]
    
    # Add to table
    row_data <- data.frame(
      Theta = theta,
      P_Cat0 = round(cat_probs[1], 3),
      P_Cat1 = round(cat_probs[2], 3),
      P_Cat2 = round(cat_probs[3], 3),
      P_Cat3 = round(cat_probs[4], 3),
      P_Cat4 = round(cat_probs[5], 3)
    )
    
    probability_table <- rbind(probability_table, row_data)
  }
  
  cat("\nCategory Probabilities by Ability Level:\n")
  print(probability_table)
  
  # Verify probabilities sum to 1
  row_sums <- rowSums(probability_table[, -1])
  cat("\nRow sums (should all be 1.000):", paste(round(row_sums, 3), collapse = ", "), "\n")
  
  # Expected score calculation
  cat("\n=== EXPECTED SCORE FUNCTION ===\n")
  
  expected_scores <- numeric(length(theta_examples))
  
  for(i in 1:length(theta_examples)) {
    theta <- theta_examples[i]
    cat_probs <- grm_probability(theta, a_example, b_example)
    expected_scores[i] <- sum(cat_probs * (0:(grm_data$n_categories - 1)))
  }
  
  expected_score_table <- data.frame(
    Theta = theta_examples,
    Expected_Score = round(expected_scores, 3),
    Min_Possible = 0,
    Max_Possible = grm_data$n_categories - 1
  )
  
  cat("Expected Score Function:\n")
  print(expected_score_table)
  
  cat("\nExpected Score Properties:\n")
  cat("• Monotonically increasing with θ\n")
  cat("• Bounded between 0 and", grm_data$n_categories - 1, "\n")
  cat("• Smooth, continuous function\n")
  cat("• Steepness determined by discrimination parameter\n")
  
  return(list(
    grm_data = grm_data,
    example_item = example_item,
    example_params = list(discrimination = a_example, thresholds = b_example),
    probability_table = probability_table,
    expected_score_table = expected_score_table
  ))
}

# Run GRM mathematical demonstration
grm_math_results <- demonstrate_grm_mathematics()
## === GRM MATHEMATICAL DEMONSTRASIYASI ===
## GRM Example Dataset:
## Items: 6 
## Persons: 1500 
## Categories per item: 5 (0-4)
## Response range: 0 4 
## 
## === GRM MATHEMATICAL FORMULATION ===
## **CUMULATIVE PROBABILITY MODEL:** 
##  
## P*ik(θ) = P(Xij ≥ k | θi) = exp[ai(θi - bik)] / [1 + exp[ai(θi - bik)]] 
##  
## Where: 
## • P*ik(θ) = Cumulative probability of responding in category k or higher 
## • ai = Discrimination parameter for item i 
## • bik = Threshold parameter k for item i 
## • θi = Ability level of person i 
##  
## **CATEGORY PROBABILITIES:** 
##  
## P(Xij = k | θi) = P*ik(θ) - P*i(k+1)(θ) 
##  
## Special cases: 
## • P(Xij = 0 | θi) = 1 - P*i1(θ) 
## • P(Xij = mi | θi) = P*imi(θ)  [highest category] 
##  
## **CONSTRAINTS:** 
##  
## • bi1 < bi2 < ... < bi(mi) [Ordered thresholds] 
## • Σk P(Xij = k | θi) = 1 [Probabilities sum to 1] 
## • All P(Xij = k | θi) ≥ 0 [Non-negative probabilities] 
## 
## === NUMERICAL EXAMPLE ===
## Example Item 1 :
## Discrimination (a): 0.811 
## Thresholds (b): -0.935, -0.217, 0.068, 0.742 
## 
## Category Probabilities by Ability Level:
##   Theta P_Cat0 P_Cat1 P_Cat2 P_Cat3 P_Cat4
## 1    -2  0.704  0.106  0.033  0.060  0.098
## 2    -1  0.513  0.141  0.050  0.100  0.196
## 3     0  0.319  0.137  0.057  0.132  0.354
## 4     1  0.172  0.099  0.048  0.128  0.552
## 5     2  0.085  0.057  0.030  0.092  0.735
## 
## Row sums (should all be 1.000): 1.001, 1, 0.999, 0.999, 0.999 
## 
## === EXPECTED SCORE FUNCTION ===
## Expected Score Function:
##   Theta Expected_Score Min_Possible Max_Possible
## 1    -2          0.742            0            4
## 2    -1          1.325            0            4
## 3     0          2.065            0            4
## 4     1          2.789            0            4
## 5     2          3.336            0            4
## 
## Expected Score Properties:
## • Monotonically increasing with θ
## • Bounded between 0 and 4 
## • Smooth, continuous function
## • Steepness determined by discrimination parameter
cat("\n=== GRM MODEL IDENTIFICATION ===\n")
## 
## === GRM MODEL IDENTIFICATION ===
# Model identification discussion
identification_constraints <- c(
  "**IDENTIFICATION CONSTRAINTS IN GRM:**",
  "",
  "1. **Ability Scale:**",
  "   • Mean(θ) = 0, Var(θ) = 1 [Standard normal distribution]",
  "   • Or fix two anchor items",
  "",
  "2. **Threshold Ordering:**",
  "   • bi1 < bi2 < ... < bimi for each item i",
  "   • Ensures proper category ordering",
  "",
  "3. **Parameter Constraints:**",
  "   • Discrimination parameters ai > 0",
  "   • No upper bound on discrimination",
  "   • Threshold parameters unbounded but ordered",
  "",
  "**COMMON PARAMETERIZATIONS:**",
  "",
  "1. **Slope-Intercept (SI):**",
  "   • P*ik = exp(ai*θ + cik) / [1 + exp(ai*θ + cik)]",
  "   • Where cik = -ai*bik",
  "",
  "2. **Slope-Threshold (ST):**", 
  "   • P*ik = exp[ai(θ - bik)] / [1 + exp[ai(θ - bik)]]",
  "   • Most common in IRT software",
  "",
  "**REPARAMETERIZATION RELATIONSHIPS:**",
  "• SI to ST: bik = -cik/ai",
  "• ST to SI: cik = -ai*bik"
)

for(constraint in identification_constraints) {
  cat(constraint, "\n")
}
## **IDENTIFICATION CONSTRAINTS IN GRM:** 
##  
## 1. **Ability Scale:** 
##    • Mean(θ) = 0, Var(θ) = 1 [Standard normal distribution] 
##    • Or fix two anchor items 
##  
## 2. **Threshold Ordering:** 
##    • bi1 < bi2 < ... < bimi for each item i 
##    • Ensures proper category ordering 
##  
## 3. **Parameter Constraints:** 
##    • Discrimination parameters ai > 0 
##    • No upper bound on discrimination 
##    • Threshold parameters unbounded but ordered 
##  
## **COMMON PARAMETERIZATIONS:** 
##  
## 1. **Slope-Intercept (SI):** 
##    • P*ik = exp(ai*θ + cik) / [1 + exp(ai*θ + cik)] 
##    • Where cik = -ai*bik 
##  
## 2. **Slope-Threshold (ST):** 
##    • P*ik = exp[ai(θ - bik)] / [1 + exp[ai(θ - bik)]] 
##    • Most common in IRT software 
##  
## **REPARAMETERIZATION RELATIONSHIPS:** 
## • SI to ST: bik = -cik/ai 
## • ST to SI: cik = -ai*bik
cat("\n=== GRM PARAMETER INTERPRETATION ===\n")
## 
## === GRM PARAMETER INTERPRETATION ===
# Parameter interpretation guide
parameter_interpretation <- data.frame(
  Parameter = c("ai (Discrimination)", "bik (Threshold k)", "Expected Score", "Information"),
  
  Interpretation = c("Item's ability to differentiate between ability levels",
                    "Ability level where P(X≥k) = 0.5",
                    "Average response for given ability",
                    "Precision of ability estimation"),
  
  Range = c("0 < ai < ∞ (typically 0.5-3.0)",
           "-∞ < bik < ∞ (typically -3 to +3)",
           "0 ≤ E[X|θ] ≤ m (max category)",
           "0 ≤ I(θ) < ∞"),
  
  High_Values = c("Steep expected score function",
                 "High ability needed for category k+",
                 "High ability level",
                 "Precise ability estimation"),
  
  Low_Values = c("Flat expected score function",
                "Low ability sufficient for category k+",
                "Low ability level",
                "Imprecise ability estimation")
)

print(parameter_interpretation)
##             Parameter                                         Interpretation
## 1 ai (Discrimination) Item's ability to differentiate between ability levels
## 2   bik (Threshold k)                       Ability level where P(X≥k) = 0.5
## 3      Expected Score                     Average response for given ability
## 4         Information                        Precision of ability estimation
##                               Range                         High_Values
## 1    0 < ai < ∞ (typically 0.5-3.0)       Steep expected score function
## 2 -∞ < bik < ∞ (typically -3 to +3) High ability needed for category k+
## 3     0 ≤ E[X|θ] ≤ m (max category)                  High ability level
## 4                      0 ≤ I(θ) < ∞          Precise ability estimation
##                               Low_Values
## 1           Flat expected score function
## 2 Low ability sufficient for category k+
## 3                      Low ability level
## 4           Imprecise ability estimation

3 Category Response Functions

3.1 Understanding CRFs və Visualization

cat("=== CATEGORY RESPONSE FUNCTIONS ===\n\n")
## === CATEGORY RESPONSE FUNCTIONS ===
cat("Category Response Functions (CRFs): Hər kateqoriya üçün ehtimal əyriləri\n")
## Category Response Functions (CRFs): Hər kateqoriya üçün ehtimal əyriləri
cat("• Müxtəlif ability səviyyələrində kateqoriya seçim ehtimalları\n")
## • Müxtəlif ability səviyyələrində kateqoriya seçim ehtimalları
cat("• Item charakteristikasının qrafik təmsili\n")
## • Item charakteristikasının qrafik təmsili
cat("• Threshold parameters və discrimination təsiri\n")
## • Threshold parameters və discrimination təsiri
cat("• Model fit və item functioning qiymətləndirməsi\n\n")
## • Model fit və item functioning qiymətləndirməsi
# Comprehensive CRF demonstration
demonstrate_category_response_functions <- function() {
  
  cat("=== CATEGORY RESPONSE FUNCTIONS DEMONSTRASIYASI ===\n")
  
  # Create examples with different parameter configurations
  theta_range <- seq(-4, 4, 0.1)
  
  # Example 1: Well-functioning item
  example1_params <- list(
    discrimination = 1.5,
    thresholds = c(-1.5, -0.5, 0.5, 1.5),
    title = "Well-Functioning Item"
  )
  
  # Example 2: Low discrimination item
  example2_params <- list(
    discrimination = 0.8,
    thresholds = c(-1.5, -0.5, 0.5, 1.5),
    title = "Low Discrimination Item"
  )
  
  # Example 3: Highly discriminating item
  example3_params <- list(
    discrimination = 2.5,
    thresholds = c(-1.5, -0.5, 0.5, 1.5),
    title = "High Discrimination Item"
  )
  
  # Example 4: Poorly spaced thresholds
  example4_params <- list(
    discrimination = 1.5,
    thresholds = c(-2.0, -1.8, -1.6, -1.4),
    title = "Poorly Spaced Thresholds"
  )
  
  examples <- list(example1_params, example2_params, example3_params, example4_params)
  
  # Calculate CRFs for each example
  calculate_crfs <- function(params, theta_range) {
    
    crf_data <- data.frame()
    
    for(theta in theta_range) {
      
      cat_probs <- grm_probability(theta, params$discrimination, params$thresholds)
      
      for(cat in 1:length(cat_probs)) {
        crf_data <- rbind(crf_data, data.frame(
          Theta = theta,
          Category = paste0("Cat", cat - 1),
          Probability = cat_probs[cat],
          Item = params$title
        ))
      }
    }
    
    return(crf_data)
  }
  
  # Generate CRF data for all examples
  all_crf_data <- data.frame()
  
  for(i in 1:length(examples)) {
    example_crf <- calculate_crfs(examples[[i]], theta_range)
    all_crf_data <- rbind(all_crf_data, example_crf)
  }
  
  # Create CRF plots
  if(require(ggplot2, quietly = TRUE)) {
    
    cat("Creating Category Response Function plots...\n")
    
    # Plot for each example
    for(i in 1:length(examples)) {
      
      example_data <- all_crf_data[all_crf_data$Item == examples[[i]]$title, ]
      
      p <- ggplot(example_data, aes(x = Theta, y = Probability, color = Category)) +
        geom_line(size = 1.2) +
        labs(title = paste("Category Response Functions:", examples[[i]]$title),
             subtitle = paste("Discrimination =", examples[[i]]$discrimination, 
                              "| Thresholds =", paste(examples[[i]]$thresholds, collapse = ", ")),
             x = "Ability (θ)",
             y = "Probability",
             color = "Category") +
        scale_color_brewer(type = "qual", palette = "Set1") +
        theme_minimal() +
        theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
              plot.subtitle = element_text(hjust = 0.5, size = 11),
              legend.position = "bottom") +
        ylim(0, 1) +
        xlim(-4, 4)
      
      print(p)
    }
  }
  
  # Analyze CRF characteristics
  cat("\n=== CRF CHARACTERISTICS ANALYSIS ===\n")
  
  analyze_crf_properties <- function(params, theta_range) {
    
    # Calculate modal categories across theta range
    modal_categories <- numeric(length(theta_range))
    max_probabilities <- numeric(length(theta_range))
    
    for(i in 1:length(theta_range)) {
      cat_probs <- grm_probability(theta_range[i], params$discrimination, params$thresholds)
      modal_categories[i] <- which.max(cat_probs) - 1  # 0-indexed
      max_probabilities[i] <- max(cat_probs)
    }
    
    # Find threshold crossing points (where adjacent categories have equal probability)
    crossing_points <- numeric(length(params$thresholds))
    
    for(k in 1:length(params$thresholds)) {
      # Find theta where P(X=k-1) = P(X=k)
      theta_k <- params$thresholds[k]  # Approximate crossing point
      
      # Refine with optimization
      objective <- function(theta) {
        cat_probs <- grm_probability(theta, params$discrimination, params$thresholds)
        abs(cat_probs[k] - cat_probs[k+1])
      }
      
      tryCatch({
        result <- optimize(objective, interval = c(theta_k - 1, theta_k + 1))
        crossing_points[k] <- result$minimum
      }, error = function(e) {
        crossing_points[k] <- theta_k
      })
    }
    
    # Category usage statistics
    category_usage <- table(modal_categories)
    n_categories_used <- length(category_usage)
    
    properties <- list(
      modal_categories = modal_categories,
      max_probabilities = max_probabilities,
      crossing_points = crossing_points,
      category_usage = category_usage,
      n_categories_used = n_categories_used,
      min_max_probability = min(max_probabilities),
      effective_range = range(theta_range[max_probabilities > 0.4])
    )
    
    return(properties)
  }
  
  # Analyze each example
  for(i in 1:length(examples)) {
    
    cat("Analysis for", examples[[i]]$title, ":\n")
    
    properties <- analyze_crf_properties(examples[[i]], theta_range)
    
    cat("Categories used:", paste(names(properties$category_usage), collapse = ", "), "\n")
    cat("Crossing points:", paste(round(properties$crossing_points, 2), collapse = ", "), "\n")
    cat("Min modal probability:", round(properties$min_max_probability, 3), "\n")
    cat("Effective range:", paste(round(properties$effective_range, 2), collapse = " to "), "\n")
    
    # Quality assessment
    quality_score <- 0
    
    # All categories used
    if(properties$n_categories_used == length(examples[[i]]$thresholds) + 1) {
      quality_score <- quality_score + 2
    }
    
    # High discrimination
    if(examples[[i]]$discrimination > 1.0) {
      quality_score <- quality_score + 1
    }
    
    # Good modal probabilities
    if(properties$min_max_probability > 0.3) {
      quality_score <- quality_score + 1
    }
    
    # Well-spaced thresholds
    threshold_spacings <- diff(examples[[i]]$thresholds)
    if(all(threshold_spacings > 0.5) && all(threshold_spacings < 2.5)) {
      quality_score <- quality_score + 1
    }
    
    quality_assessment <- ifelse(quality_score >= 4, "Excellent",
                                ifelse(quality_score >= 3, "Good",
                                      ifelse(quality_score >= 2, "Fair", "Poor")))
    
    cat("Quality assessment:", quality_assessment, "(", quality_score, "/5 )\n\n")
  }
  
  return(list(
    examples = examples,
    crf_data = all_crf_data,
    theta_range = theta_range
  ))
}

# Run CRF demonstration
crf_results <- demonstrate_category_response_functions()
## === CATEGORY RESPONSE FUNCTIONS DEMONSTRASIYASI ===
## Creating Category Response Function plots...

## 
## === CRF CHARACTERISTICS ANALYSIS ===
## Analysis for Well-Functioning Item :
## Categories used: 0, 1, 2, 3, 4 
## Crossing points: -1.11, -0.5, 0.5, 1.11 
## Min modal probability: 0.318 
## Effective range: -4 to 4 
## Quality assessment: Excellent ( 5 /5 )
## 
## Analysis for Low Discrimination Item :
## Categories used: 0, 4 
## Crossing points: -0.5, -0.5, 0.5, 0.5 
## Min modal probability: 0.231 
## Effective range: -4 to 4 
## Quality assessment: Poor ( 1 /5 )
## 
## Analysis for High Discrimination Item :
## Categories used: 0, 1, 2, 3, 4 
## Crossing points: -1.43, -0.5, 0.5, 1.43 
## Min modal probability: 0.424 
## Effective range: -4 to 4 
## Quality assessment: Excellent ( 5 /5 )
## 
## Analysis for Poorly Spaced Thresholds :
## Categories used: 0, 4 
## Crossing points: -1, -1.8, -1.6, -2.4 
## Min modal probability: 0.389 
## Effective range: -4 to 4 
## Quality assessment: Fair ( 2 /5 )
cat("\n=== CRF INTERPRETATION GUIDELINES ===\n")
## 
## === CRF INTERPRETATION GUIDELINES ===
# CRF interpretation guidelines
crf_guidelines <- c(
  "**IDEAL CRF CHARACTERISTICS:**",
  "",
  "1. **Category Usage:**",
  "   • All categories should be modal somewhere on θ scale",
  "   • Each category should have peak probability > 0.3",
  "   • No 'dead' categories that are never most likely",
  "",
  "2. **Threshold Spacing:**",
  "   • Adjacent thresholds spaced 0.5 to 2.0 logits apart",
  "   • Equal spacing often optimal for rating scales",
  "   • Avoid clustering of thresholds",
  "",
  "3. **Discrimination Effects:**",
  "   • Higher discrimination → steeper CRFs",
  "   • Sharper transitions between categories",
  "   • More precise ability discrimination",
  "",
  "4. **Troubleshooting Poor CRFs:**",
  "   • Flat CRFs → increase discrimination or collapse categories",
  "   • Unused categories → adjust thresholds or reduce categories",
  "   • Disordered thresholds → violates model assumptions"
)

for(guideline in crf_guidelines) {
  cat(guideline, "\n")
}
## **IDEAL CRF CHARACTERISTICS:** 
##  
## 1. **Category Usage:** 
##    • All categories should be modal somewhere on θ scale 
##    • Each category should have peak probability > 0.3 
##    • No 'dead' categories that are never most likely 
##  
## 2. **Threshold Spacing:** 
##    • Adjacent thresholds spaced 0.5 to 2.0 logits apart 
##    • Equal spacing often optimal for rating scales 
##    • Avoid clustering of thresholds 
##  
## 3. **Discrimination Effects:** 
##    • Higher discrimination → steeper CRFs 
##    • Sharper transitions between categories 
##    • More precise ability discrimination 
##  
## 4. **Troubleshooting Poor CRFs:** 
##    • Flat CRFs → increase discrimination or collapse categories 
##    • Unused categories → adjust thresholds or reduce categories 
##    • Disordered thresholds → violates model assumptions
cat("\n=== CRF QUALITY INDICATORS ===\n")
## 
## === CRF QUALITY INDICATORS ===
# Quality indicators table
quality_indicators <- data.frame(
  Indicator = c("Category Usage", "Peak Probabilities", "Threshold Spacing", 
                "Discrimination Level", "Effective Range", "Model Fit"),
  
  Good_Range = c("All categories used", "> 0.30", "0.5 - 2.0 logits", 
                 "> 1.0", "Wide coverage", "Adequate"),
  
  Warning_Signs = c("Unused categories", "< 0.20", "< 0.5 or > 2.5", 
                   "< 0.7", "Narrow range", "Poor fit"),
  
  Remedial_Actions = c("Collapse categories", "Adjust thresholds", "Redistribute thresholds", 
                      "Check item quality", "Review content", "Model modification"),
  
  Software_Checks = c("Modal category plots", "CRF maximum values", "Threshold differences", 
                     "Parameter estimates", "Information plots", "Fit statistics")
)

print(quality_indicators)
##              Indicator          Good_Range     Warning_Signs
## 1       Category Usage All categories used Unused categories
## 2   Peak Probabilities              > 0.30            < 0.20
## 3    Threshold Spacing    0.5 - 2.0 logits    < 0.5 or > 2.5
## 4 Discrimination Level               > 1.0             < 0.7
## 5      Effective Range       Wide coverage      Narrow range
## 6            Model Fit            Adequate          Poor fit
##          Remedial_Actions       Software_Checks
## 1     Collapse categories  Modal category plots
## 2       Adjust thresholds    CRF maximum values
## 3 Redistribute thresholds Threshold differences
## 4      Check item quality   Parameter estimates
## 5          Review content     Information plots
## 6      Model modification        Fit statistics

4 Threshold Parameters

4.1 Understanding və Interpretation

cat("=== THRESHOLD PARAMETERS ===\n\n")
## === THRESHOLD PARAMETERS ===
cat("Threshold Parameters: GRM-də kateqoriya keçidlərini müəyyən edən parameterlər\n")
## Threshold Parameters: GRM-də kateqoriya keçidlərini müəyyən edən parameterlər
cat("• bik: i item-inin k threshold-u\n")
## • bik: i item-inin k threshold-u
cat("• P(X ≥ k | θ = bik) = 0.5 interpretasiyası\n")
## • P(X ≥ k | θ = bik) = 0.5 interpretasiyası
cat("• Threshold ordering və constraints\n")
## • Threshold ordering və constraints
cat("• Praktik interpretasiya və meaningful differences\n\n")
## • Praktik interpretasiya və meaningful differences
# Comprehensive threshold parameters demonstration
demonstrate_threshold_parameters <- function() {
  
  cat("=== THRESHOLD PARAMETERS DEMONSTRASIYASI ===\n")
  
  # Create examples with different threshold configurations
  theta_range <- seq(-3, 3, 0.1)
  
  # Example scenarios for threshold analysis
  threshold_scenarios <- list(
    
    "Well_Spaced" = list(
      discrimination = 1.5,
      thresholds = c(-1.5, -0.5, 0.5, 1.5),
      description = "Well-spaced thresholds"
    ),
    
    "Close_Thresholds" = list(
      discrimination = 1.5,
      thresholds = c(-0.3, -0.1, 0.1, 0.3),
      description = "Closely spaced thresholds"
    ),
    
    "Wide_Thresholds" = list(
      discrimination = 1.5,
      thresholds = c(-2.5, -1.0, 1.0, 2.5),
      description = "Widely spaced thresholds"
    ),
    
    "Unequal_Spacing" = list(
      discrimination = 1.5,
      thresholds = c(-2.0, -1.8, 0.5, 2.0),
      description = "Unequally spaced thresholds"
    )
  )
  
  cat("Threshold Scenarios Created:\n")
  for(name in names(threshold_scenarios)) {
    scenario <- threshold_scenarios[[name]]
    cat(name, ":", scenario$description, "\n")
    cat("  Thresholds:", paste(scenario$thresholds, collapse = ", "), "\n")
  }
  cat("\n")
  
  # Analyze threshold properties
  analyze_threshold_properties <- function(scenario_name, params) {
    
    cat("=== ANALYSIS:", scenario_name, "===\n")
    
    # Basic threshold statistics
    thresholds <- params$thresholds
    n_thresholds <- length(thresholds)
    
    # Threshold spacing
    threshold_spacing <- diff(thresholds)
    
    cat("Threshold Analysis:\n")
    cat("Number of thresholds:", n_thresholds, "\n")
    cat("Threshold values:", paste(round(thresholds, 2), collapse = ", "), "\n")
    cat("Threshold spacings:", paste(round(threshold_spacing, 2), collapse = ", "), "\n")
    cat("Mean spacing:", round(mean(threshold_spacing), 2), "\n")
    cat("SD spacing:", round(sd(threshold_spacing), 2), "\n")
    
    # Expected score at thresholds
    expected_scores_at_thresholds <- numeric(n_thresholds)
    
    for(i in 1:n_thresholds) {
      theta_threshold <- thresholds[i]
      cat_probs <- grm_probability(theta_threshold, params$discrimination, thresholds)
      expected_scores_at_thresholds[i] <- sum(cat_probs * (0:(length(cat_probs) - 1)))
    }
    
    cat("Expected scores at thresholds:", paste(round(expected_scores_at_thresholds, 2), collapse = ", "), "\n")
    
    # Calculate cumulative probabilities at thresholds
    cumulative_probs_at_thresholds <- matrix(NA, n_thresholds, n_thresholds)
    
    for(i in 1:n_thresholds) {
      theta_threshold <- thresholds[i]
      
      for(k in 1:n_thresholds) {
        cum_prob <- exp(params$discrimination * (theta_threshold - thresholds[k])) / 
                    (1 + exp(params$discrimination * (theta_threshold - thresholds[k])))
        cumulative_probs_at_thresholds[i, k] <- cum_prob
      }
    }
    
    cat("Cumulative probabilities at own threshold (should be ≈ 0.5):\n")
    for(i in 1:n_thresholds) {
      cat("  Threshold", i, ":", round(cumulative_probs_at_thresholds[i, i], 3), "\n")
    }
    
    # Category probabilities at each threshold
    cat("\nCategory probabilities at each threshold:\n")
    
    for(i in 1:n_thresholds) {
      theta_threshold <- thresholds[i]
      cat_probs <- grm_probability(theta_threshold, params$discrimination, thresholds)
      
      cat("At threshold", i, "(θ =", round(theta_threshold, 2), "):\n")
      for(k in 1:length(cat_probs)) {
        cat("  P(X =", k-1, ") =", round(cat_probs[k], 3), "\n")
      }
    }
    
    # Quality assessment
    quality_indicators <- list()
    
    # Spacing uniformity
    spacing_cv <- sd(threshold_spacing) / mean(threshold_spacing)
    quality_indicators$spacing_uniformity <- ifelse(spacing_cv < 0.3, "Good", "Poor")
    
    # Adequate spacing
    min_spacing <- min(threshold_spacing)
    quality_indicators$adequate_spacing <- ifelse(min_spacing > 0.5, "Good", "Concern")
    
    # Not too wide
    max_spacing <- max(threshold_spacing)
    quality_indicators$not_too_wide <- ifelse(max_spacing < 2.5, "Good", "Concern")
    
    # Expected score progression
    score_increases <- diff(expected_scores_at_thresholds)
    quality_indicators$score_progression <- ifelse(all(score_increases > 0), "Good", "Problem")
    
    cat("\nQuality Assessment:\n")
    for(indicator in names(quality_indicators)) {
      cat(indicator, ":", quality_indicators[[indicator]], "\n")
    }
    
    cat("\n")
    
    return(list(
      thresholds = thresholds,
      spacing = threshold_spacing,
      expected_scores = expected_scores_at_thresholds,
      quality = quality_indicators
    ))
  }
  
  # Analyze all scenarios
  scenario_analyses <- list()
  
  for(scenario_name in names(threshold_scenarios)) {
    analysis <- analyze_threshold_properties(scenario_name, threshold_scenarios[[scenario_name]])
    scenario_analyses[[scenario_name]] <- analysis
  }
  
  # Create threshold visualization
  if(require(ggplot2, quietly = TRUE)) {
    
    cat("=== THRESHOLD VISUALIZATION ===\n")
    
    # Create cumulative probability plots
    for(scenario_name in names(threshold_scenarios)) {
      
      params <- threshold_scenarios[[scenario_name]]
      
      # Calculate cumulative probabilities
      cum_prob_data <- data.frame()
      
      for(theta in theta_range) {
        for(k in 1:length(params$thresholds)) {
          cum_prob <- exp(params$discrimination * (theta - params$thresholds[k])) / 
                      (1 + exp(params$discrimination * (theta - params$thresholds[k])))
          
          cum_prob_data <- rbind(cum_prob_data, data.frame(
            Theta = theta,
            Threshold = paste0("P*(≥", k, ")"),
            Cumulative_Probability = cum_prob,
            Scenario = params$description
          ))
        }
      }
      
      # Create plot
      p <- ggplot(cum_prob_data, aes(x = Theta, y = Cumulative_Probability, color = Threshold)) +
        geom_line(size = 1.2) +
        geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.5) +
        geom_vline(xintercept = params$thresholds, linetype = "dotted", alpha = 0.7) +
        labs(title = paste("Cumulative Probabilities:", params$description),
             subtitle = paste("Thresholds:", paste(round(params$thresholds, 2), collapse = ", ")),
             x = "Ability (θ)",
             y = "Cumulative Probability P*(X ≥ k)",
             color = "Threshold") +
        scale_color_brewer(type = "qual", palette = "Dark2") +
        theme_minimal() +
        theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
              plot.subtitle = element_text(hjust = 0.5, size = 11),
              legend.position = "bottom") +
        ylim(0, 1) +
        xlim(-3, 3)
      
      print(p)
    }
  }
  
  # Threshold interpretation guidelines
  cat("=== THRESHOLD INTERPRETATION GUIDELINES ===\n")
  
  interpretation_guidelines <- c(
    "**THRESHOLD PARAMETER MEANING:**",
    "",
    "bik = Ability level where P(X ≥ k | θ) = 0.5",
    "",
    "• Threshold 1 (bi1): 50% chance of scoring 1 or higher",
    "• Threshold 2 (bi2): 50% chance of scoring 2 or higher", 
    "• And so on...",
    "",
    "**PRACTICAL INTERPRETATION:**",
    "",
    "• Lower thresholds → easier to achieve higher categories",
    "• Higher thresholds → harder to achieve higher categories",
    "• Threshold spacing affects category usage patterns",
    "",
    "**OPTIMAL THRESHOLD SPACING:**",
    "",
    "• Generally 0.5 to 2.0 logits apart",
    "• Equal spacing often works well for rating scales",
    "• Closer spacing → more categories used near center",
    "• Wider spacing → broader category usage",
    "",
    "**THRESHOLD PROBLEMS:**",
    "",
    "• Disordered thresholds (bi1 > bi2): Model violation",
    "• Very close thresholds (< 0.3): Categories hard to distinguish",
    "• Very distant thresholds (> 3.0): Missing intermediate responses"
  )
  
  for(guideline in interpretation_guidelines) {
    cat(guideline, "\n")
  }
  
  return(list(
    scenarios = threshold_scenarios,
    analyses = scenario_analyses
  ))
}

# Run threshold parameters demonstration
threshold_results <- demonstrate_threshold_parameters()
## === THRESHOLD PARAMETERS DEMONSTRASIYASI ===
## Threshold Scenarios Created:
## Well_Spaced : Well-spaced thresholds 
##   Thresholds: -1.5, -0.5, 0.5, 1.5 
## Close_Thresholds : Closely spaced thresholds 
##   Thresholds: -0.3, -0.1, 0.1, 0.3 
## Wide_Thresholds : Widely spaced thresholds 
##   Thresholds: -2.5, -1, 1, 2.5 
## Unequal_Spacing : Unequally spaced thresholds 
##   Thresholds: -2, -1.8, 0.5, 2 
## 
## === ANALYSIS: Well_Spaced ===
## Threshold Analysis:
## Number of thresholds: 4 
## Threshold values: -1.5, -0.5, 0.5, 1.5 
## Threshold spacings: 1, 1, 1 
## Mean spacing: 1 
## SD spacing: 0 
## Expected scores at thresholds: 0.74, 1.55, 2.45, 3.26 
## Cumulative probabilities at own threshold (should be ≈ 0.5):
##   Threshold 1 : 0.5 
##   Threshold 2 : 0.5 
##   Threshold 3 : 0.5 
##   Threshold 4 : 0.5 
## 
## Category probabilities at each threshold:
## At threshold 1 (θ = -1.5 ):
##   P(X = 0 ) = 0.5 
##   P(X = 1 ) = 0.318 
##   P(X = 2 ) = 0.135 
##   P(X = 3 ) = 0.036 
##   P(X = 4 ) = 0.011 
## At threshold 2 (θ = -0.5 ):
##   P(X = 0 ) = 0.182 
##   P(X = 1 ) = 0.318 
##   P(X = 2 ) = 0.318 
##   P(X = 3 ) = 0.135 
##   P(X = 4 ) = 0.047 
## At threshold 3 (θ = 0.5 ):
##   P(X = 0 ) = 0.047 
##   P(X = 1 ) = 0.135 
##   P(X = 2 ) = 0.318 
##   P(X = 3 ) = 0.318 
##   P(X = 4 ) = 0.182 
## At threshold 4 (θ = 1.5 ):
##   P(X = 0 ) = 0.011 
##   P(X = 1 ) = 0.036 
##   P(X = 2 ) = 0.135 
##   P(X = 3 ) = 0.318 
##   P(X = 4 ) = 0.5 
## 
## Quality Assessment:
## spacing_uniformity : Good 
## adequate_spacing : Good 
## not_too_wide : Good 
## score_progression : Good 
## 
## === ANALYSIS: Close_Thresholds ===
## Threshold Analysis:
## Number of thresholds: 4 
## Threshold values: -0.3, -0.1, 0.1, 0.3 
## Threshold spacings: 0.2, 0.2, 0.2 
## Mean spacing: 0.2 
## SD spacing: 0 
## Expected scores at thresholds: 1.57, 1.85, 2.15, 2.43 
## Cumulative probabilities at own threshold (should be ≈ 0.5):
##   Threshold 1 : 0.5 
##   Threshold 2 : 0.5 
##   Threshold 3 : 0.5 
##   Threshold 4 : 0.5 
## 
## Category probabilities at each threshold:
## At threshold 1 (θ = -0.3 ):
##   P(X = 0 ) = 0.5 
##   P(X = 1 ) = 0.074 
##   P(X = 2 ) = 0.071 
##   P(X = 3 ) = 0.065 
##   P(X = 4 ) = 0.289 
## At threshold 2 (θ = -0.1 ):
##   P(X = 0 ) = 0.426 
##   P(X = 1 ) = 0.074 
##   P(X = 2 ) = 0.074 
##   P(X = 3 ) = 0.071 
##   P(X = 4 ) = 0.354 
## At threshold 3 (θ = 0.1 ):
##   P(X = 0 ) = 0.354 
##   P(X = 1 ) = 0.071 
##   P(X = 2 ) = 0.074 
##   P(X = 3 ) = 0.074 
##   P(X = 4 ) = 0.426 
## At threshold 4 (θ = 0.3 ):
##   P(X = 0 ) = 0.289 
##   P(X = 1 ) = 0.065 
##   P(X = 2 ) = 0.071 
##   P(X = 3 ) = 0.074 
##   P(X = 4 ) = 0.5 
## 
## Quality Assessment:
## spacing_uniformity : Good 
## adequate_spacing : Concern 
## not_too_wide : Good 
## score_progression : Good 
## 
## === ANALYSIS: Wide_Thresholds ===
## Threshold Analysis:
## Number of thresholds: 4 
## Threshold values: -2.5, -1, 1, 2.5 
## Threshold spacings: 1.5, 2, 1.5 
## Mean spacing: 1.67 
## SD spacing: 0.29 
## Expected scores at thresholds: 0.6, 1.46, 2.54, 3.4 
## Cumulative probabilities at own threshold (should be ≈ 0.5):
##   Threshold 1 : 0.5 
##   Threshold 2 : 0.5 
##   Threshold 3 : 0.5 
##   Threshold 4 : 0.5 
## 
## Category probabilities at each threshold:
## At threshold 1 (θ = -2.5 ):
##   P(X = 0 ) = 0.5 
##   P(X = 1 ) = 0.405 
##   P(X = 2 ) = 0.09 
##   P(X = 3 ) = 0.005 
##   P(X = 4 ) = 0.001 
## At threshold 2 (θ = -1 ):
##   P(X = 0 ) = 0.095 
##   P(X = 1 ) = 0.405 
##   P(X = 2 ) = 0.453 
##   P(X = 3 ) = 0.042 
##   P(X = 4 ) = 0.005 
## At threshold 3 (θ = 1 ):
##   P(X = 0 ) = 0.005 
##   P(X = 1 ) = 0.042 
##   P(X = 2 ) = 0.453 
##   P(X = 3 ) = 0.405 
##   P(X = 4 ) = 0.095 
## At threshold 4 (θ = 2.5 ):
##   P(X = 0 ) = 0.001 
##   P(X = 1 ) = 0.005 
##   P(X = 2 ) = 0.09 
##   P(X = 3 ) = 0.405 
##   P(X = 4 ) = 0.5 
## 
## Quality Assessment:
## spacing_uniformity : Good 
## adequate_spacing : Good 
## not_too_wide : Good 
## score_progression : Good 
## 
## === ANALYSIS: Unequal_Spacing ===
## Threshold Analysis:
## Number of thresholds: 4 
## Threshold values: -2, -1.8, 0.5, 2 
## Threshold spacings: 0.2, 2.3, 1.5 
## Mean spacing: 1.33 
## SD spacing: 1.06 
## Expected scores at thresholds: 0.95, 1.11, 2.54, 3.4 
## Cumulative probabilities at own threshold (should be ≈ 0.5):
##   Threshold 1 : 0.5 
##   Threshold 2 : 0.5 
##   Threshold 3 : 0.5 
##   Threshold 4 : 0.5 
## 
## Category probabilities at each threshold:
## At threshold 1 (θ = -2 ):
##   P(X = 0 ) = 0.5 
##   P(X = 1 ) = 0.074 
##   P(X = 2 ) = 0.403 
##   P(X = 3 ) = 0.021 
##   P(X = 4 ) = 0.002 
## At threshold 2 (θ = -1.8 ):
##   P(X = 0 ) = 0.426 
##   P(X = 1 ) = 0.074 
##   P(X = 2 ) = 0.469 
##   P(X = 3 ) = 0.027 
##   P(X = 4 ) = 0.003 
## At threshold 3 (θ = 0.5 ):
##   P(X = 0 ) = 0.023 
##   P(X = 1 ) = 0.008 
##   P(X = 2 ) = 0.469 
##   P(X = 3 ) = 0.405 
##   P(X = 4 ) = 0.095 
## At threshold 4 (θ = 2 ):
##   P(X = 0 ) = 0.002 
##   P(X = 1 ) = 0.001 
##   P(X = 2 ) = 0.092 
##   P(X = 3 ) = 0.405 
##   P(X = 4 ) = 0.5 
## 
## Quality Assessment:
## spacing_uniformity : Poor 
## adequate_spacing : Concern 
## not_too_wide : Good 
## score_progression : Good 
## 
## === THRESHOLD VISUALIZATION ===

## === THRESHOLD INTERPRETATION GUIDELINES ===
## **THRESHOLD PARAMETER MEANING:** 
##  
## bik = Ability level where P(X ≥ k | θ) = 0.5 
##  
## • Threshold 1 (bi1): 50% chance of scoring 1 or higher 
## • Threshold 2 (bi2): 50% chance of scoring 2 or higher 
## • And so on... 
##  
## **PRACTICAL INTERPRETATION:** 
##  
## • Lower thresholds → easier to achieve higher categories 
## • Higher thresholds → harder to achieve higher categories 
## • Threshold spacing affects category usage patterns 
##  
## **OPTIMAL THRESHOLD SPACING:** 
##  
## • Generally 0.5 to 2.0 logits apart 
## • Equal spacing often works well for rating scales 
## • Closer spacing → more categories used near center 
## • Wider spacing → broader category usage 
##  
## **THRESHOLD PROBLEMS:** 
##  
## • Disordered thresholds (bi1 > bi2): Model violation 
## • Very close thresholds (< 0.3): Categories hard to distinguish 
## • Very distant thresholds (> 3.0): Missing intermediate responses
cat("\n=== THRESHOLD PARAMETER ESTIMATION ===\n")
## 
## === THRESHOLD PARAMETER ESTIMATION ===
# Parameter estimation considerations
estimation_considerations <- data.frame(
  Aspect = c("Sample Size", "Category Frequency", "Threshold Stability", 
             "Convergence", "Starting Values", "Constraints"),
  
  Requirement = c("50+ per category", "5+ responses per category", "Large N for precision", 
                 "Proper algorithm", "Reasonable initial values", "Ordered thresholds"),
  
  Common_Problems = c("Small samples", "Empty categories", "Boundary estimates", 
                     "Non-convergence", "Poor starting values", "Disordered estimates"),
  
  Solutions = c("Increase N or collapse", "Combine categories", "Add constraints", 
               "Change algorithm", "Use data-based starts", "Fix ordering"),
  
  Software_Notes = c("mirt, ltm handle automatically", "Automatic collapsing", "Bayesian priors help", 
                    "Multiple optimizers", "Data-driven defaults", "Built-in ordering")
)

print(estimation_considerations)
##                Aspect               Requirement      Common_Problems
## 1         Sample Size          50+ per category        Small samples
## 2  Category Frequency 5+ responses per category     Empty categories
## 3 Threshold Stability     Large N for precision   Boundary estimates
## 4         Convergence          Proper algorithm      Non-convergence
## 5     Starting Values Reasonable initial values Poor starting values
## 6         Constraints        Ordered thresholds Disordered estimates
##                Solutions                 Software_Notes
## 1 Increase N or collapse mirt, ltm handle automatically
## 2     Combine categories           Automatic collapsing
## 3        Add constraints           Bayesian priors help
## 4       Change algorithm            Multiple optimizers
## 5  Use data-based starts           Data-driven defaults
## 6           Fix ordering              Built-in ordering
cat("\n=== THRESHOLD REPORTING STANDARDS ===\n")
## 
## === THRESHOLD REPORTING STANDARDS ===
# Reporting standards
reporting_standards <- c(
  "**REQUIRED REPORTING:**",
  "",
  "1. **Parameter Estimates:**",
  "   • All threshold estimates with standard errors",
  "   • Discrimination parameters",
  "   • Model identification constraints used",
  "",
  "2. **Diagnostic Information:**",
  "   • Threshold spacing analysis",
  "   • Category usage frequencies",
  "   • Parameter stability across samples",
  "",
  "3. **Quality Indicators:**",
  "   • Model fit statistics",
  "   • Category response function plots",
  "   • Information function curves",
  "",
  "**INTERPRETATION GUIDELINES:**",
  "• Relate thresholds to practical meaning",
  "• Discuss category spacing implications",
  "• Address any problematic patterns",
  "• Provide guidance for score interpretation"
)

for(standard in reporting_standards) {
  cat(standard, "\n")
}
## **REQUIRED REPORTING:** 
##  
## 1. **Parameter Estimates:** 
##    • All threshold estimates with standard errors 
##    • Discrimination parameters 
##    • Model identification constraints used 
##  
## 2. **Diagnostic Information:** 
##    • Threshold spacing analysis 
##    • Category usage frequencies 
##    • Parameter stability across samples 
##  
## 3. **Quality Indicators:** 
##    • Model fit statistics 
##    • Category response function plots 
##    • Information function curves 
##  
## **INTERPRETATION GUIDELINES:** 
## • Relate thresholds to practical meaning 
## • Discuss category spacing implications 
## • Address any problematic patterns 
## • Provide guidance for score interpretation

5 Rating Scale Model və Partial Credit Model

5.1 Alternative Polytomous IRT Models

cat("=== RATING SCALE MODEL VƏ PARTIAL CREDIT MODEL ===\n\n")
## === RATING SCALE MODEL VƏ PARTIAL CREDIT MODEL ===
cat("Alternative Polytomous Models: GRM-ə alternativ yanaşmalar\n")
## Alternative Polytomous Models: GRM-ə alternativ yanaşmalar
cat("• Rating Scale Model (RSM): Ümumi threshold structure\n")
## • Rating Scale Model (RSM): Ümumi threshold structure
cat("• Partial Credit Model (PCM): Item-specific thresholds\n")
## • Partial Credit Model (PCM): Item-specific thresholds
cat("• Generalized Partial Credit Model (GPCM): PCM + discrimination\n")
## • Generalized Partial Credit Model (GPCM): PCM + discrimination
cat("• Model comparison və selection criteria\n\n")
## • Model comparison və selection criteria
# Comprehensive RSM and PCM demonstration
demonstrate_rsm_pcm_models <- function() {
  
  cat("=== RSM VƏ PCM MODELS DEMONSTRASIYASI ===\n")
  
  # Generate data suitable for different models
  set.seed(1801)
  n_items <- 8
  n_persons <- 1200
  n_categories <- 4  # 0, 1, 2, 3
  
  # Generate person abilities
  theta <- rnorm(n_persons, 0, 1)
  
  # Generate item difficulties (location parameters)
  item_difficulties <- rnorm(n_items, 0, 1)
  
  # Model 1: Rating Scale Model (RSM)
  # Common threshold structure across items
  common_thresholds <- c(-1.2, 0.0, 1.2)  # τ1, τ2, τ3
  
  # Model 2: Partial Credit Model (PCM)  
  # Item-specific thresholds
  item_thresholds <- matrix(NA, n_items, n_categories - 1)
  for(i in 1:n_items) {
    # Generate ordered thresholds for each item
    item_thresholds[i, ] <- sort(rnorm(n_categories - 1, 0, 0.8))
  }
  
  cat("Model Configurations:\n")
  cat("RSM - Common thresholds:", paste(round(common_thresholds, 2), collapse = ", "), "\n")
  cat("PCM - Item-specific thresholds (first 3 items):\n")
  for(i in 1:3) {
    cat("  Item", i, ":", paste(round(item_thresholds[i, ], 2), collapse = ", "), "\n")
  }
  cat("\n")
  
  # RSM probability function
  rsm_probability <- function(theta, item_difficulty, common_thresholds) {
    
    n_cats <- length(common_thresholds) + 1
    
    # Calculate numerators for each category
    numerators <- numeric(n_cats)
    
    for(k in 1:n_cats) {
      if(k == 1) {
        numerators[k] <- 1  # Category 0
      } else {
        sum_thresholds <- sum(common_thresholds[1:(k-1)])
        numerators[k] <- exp((k-1) * theta - (k-1) * item_difficulty - sum_thresholds)
      }
    }
    
    # Denominator
    denominator <- sum(numerators)
    
    # Category probabilities
    probabilities <- numerators / denominator
    
    return(probabilities)
  }
  
  # PCM probability function
  pcm_probability <- function(theta, item_difficulty, item_thresholds) {
    
    n_cats <- length(item_thresholds) + 1
    
    # Calculate numerators for each category
    numerators <- numeric(n_cats)
    
    for(k in 1:n_cats) {
      if(k == 1) {
        numerators[k] <- 1  # Category 0
      } else {
        sum_thresholds <- sum(item_thresholds[1:(k-1)])
        numerators[k] <- exp((k-1) * theta - (k-1) * item_difficulty - sum_thresholds)
      }
    }
    
    # Denominator
    denominator <- sum(numerators)
    
    # Category probabilities
    probabilities <- numerators / denominator
    
    return(probabilities)
  }
  
  # Generate data under RSM
  rsm_responses <- matrix(NA, n_persons, n_items)
  
  for(i in 1:n_persons) {
    for(j in 1:n_items) {
      cat_probs <- rsm_probability(theta[i], item_difficulties[j], common_thresholds)
      rsm_responses[i, j] <- sample(0:(n_categories-1), 1, prob = cat_probs)
    }
  }
  
  # Generate data under PCM
  pcm_responses <- matrix(NA, n_persons, n_items)
  
  for(i in 1:n_persons) {
    for(j in 1:n_items) {
      cat_probs <- pcm_probability(theta[i], item_difficulties[j], item_thresholds[j, ])
      pcm_responses[i, j] <- sample(0:(n_categories-1), 1, prob = cat_probs)
    }
  }
  
  colnames(rsm_responses) <- paste0("Item", 1:n_items)
  colnames(pcm_responses) <- paste0("Item", 1:n_items)
  
  cat("Response Data Generated:\n")
  cat("RSM data range:", range(rsm_responses), "\n")
  cat("PCM data range:", range(pcm_responses), "\n")
  
  # Descriptive analysis
  cat("\n=== DESCRIPTIVE ANALYSIS ===\n")
  
  # Category usage patterns
  analyze_category_usage <- function(responses, model_name) {
    
    cat("Category Usage -", model_name, ":\n")
    
    overall_usage <- table(factor(as.vector(responses), levels = 0:(n_categories-1)))
    cat("Overall:", paste(overall_usage, collapse = " / "), 
        "(", paste(round(prop.table(overall_usage) * 100, 1), collapse = "% / "), "%)\n")
    
    # Item-specific usage
    item_usage <- data.frame()
    
    for(j in 1:n_items) {
      item_cats <- table(factor(responses[, j], levels = 0:(n_categories-1)))
      item_usage <- rbind(item_usage, data.frame(
        Item = j,
        Cat0 = item_cats[1],
        Cat1 = item_cats[2], 
        Cat2 = item_cats[3],
        Cat3 = item_cats[4],
        Mean = round(mean(responses[, j]), 2),
        SD = round(sd(responses[, j]), 2)
      ))
    }
    
    cat("Item-level usage (first 4 items):\n")
    print(item_usage[1:4, ])
    
    return(item_usage)
  }
  
  rsm_usage <- analyze_category_usage(rsm_responses, "RSM")
  cat("\n")
  pcm_usage <- analyze_category_usage(pcm_responses, "PCM")
  
  # Model comparison using available packages
  cat("\n=== MODEL FITTING VƏ COMPARISON ===\n")
  
  if(require(eRm, quietly = TRUE)) {
    
    cat("Fitting models with eRm package...\n")
    
    tryCatch({
      # Fit RSM
      rsm_fit <- RSM(rsm_responses)
      cat("RSM fitted successfully\n")
      
      # Fit PCM  
      pcm_fit <- PCM(pcm_responses)
      cat("PCM fitted successfully\n")
      
      # Model comparison would go here
      # Note: eRm uses different parameterization than our simulation
      
    }, error = function(e) {
      cat("eRm fitting failed:", e$message, "\n")
    })
  }
  
  if(require(mirt, quietly = TRUE)) {
    
    cat("\nFitting models with mirt package...\n")
    
    tryCatch({
      # Fit RSM (constrained GPCM)
      rsm_mirt <- mirt(rsm_responses, 1, itemtype = "rsm", verbose = FALSE)
      cat("RSM fitted with mirt\n")
      
      # Fit PCM (partial credit model)
      pcm_mirt <- mirt(pcm_responses, 1, itemtype = "Rasch", verbose = FALSE)
      cat("PCM fitted with mirt\n")
      
      # Fit GPCM (generalized partial credit)
      gpcm_mirt <- mirt(pcm_responses, 1, itemtype = "gpcm", verbose = FALSE)
      cat("GPCM fitted with mirt\n")
      
      # Extract fit statistics
      rsm_fit_stats <- data.frame(
        Model = "RSM",
        AIC = extract.mirt(rsm_mirt, "AIC"),
        BIC = extract.mirt(rsm_mirt, "BIC"),
        LogLik = extract.mirt(rsm_mirt, "logLik")
      )
      
      pcm_fit_stats <- data.frame(
        Model = "PCM", 
        AIC = extract.mirt(pcm_mirt, "AIC"),
        BIC = extract.mirt(pcm_mirt, "BIC"),
        LogLik = extract.mirt(pcm_mirt, "logLik")
      )
      
      gpcm_fit_stats <- data.frame(
        Model = "GPCM",
        AIC = extract.mirt(gpcm_mirt, "AIC"),
        BIC = extract.mirt(gpcm_mirt, "BIC"),
        LogLik = extract.mirt(gpcm_mirt, "logLik")
      )
      
      fit_comparison <- rbind(rsm_fit_stats, pcm_fit_stats, gpcm_fit_stats)
      fit_comparison$Delta_AIC <- fit_comparison$AIC - min(fit_comparison$AIC)
      fit_comparison$Delta_BIC <- fit_comparison$BIC - min(fit_comparison$BIC)
      
      cat("\nModel Fit Comparison:\n")
      print(fit_comparison)
      
      # Parameter comparison
      cat("\nParameter Comparison (first 3 items):\n")
      
      rsm_coefs <- coef(rsm_mirt, simplify = TRUE)
      pcm_coefs <- coef(pcm_mirt, simplify = TRUE)
      gpcm_coefs <- coef(gpcm_mirt, simplify = TRUE)
      
      cat("RSM item parameters (sample):\n")
      print(round(rsm_coefs$items[1:3, ], 3))
      
      cat("\nGPCM item parameters (sample):\n")
      print(round(gpcm_coefs$items[1:3, ], 3))
      
    }, error = function(e) {
      cat("mirt fitting failed:", e$message, "\n")
      fit_comparison <- NULL
    })
  }
  
  # Model characteristic comparison
  cat("\n=== MODEL CHARACTERISTICS COMPARISON ===\n")
  
  model_comparison <- data.frame(
    Aspect = c("Parameterization", "Threshold Structure", "Discrimination", 
              "Complexity", "Assumptions", "Best For"),
    
    GRM = c("Cumulative logits", "Item-specific ordered", "Item-specific", 
           "Most complex", "Unidimensional", "Likert scales"),
    
    RSM = c("Adjacent categories", "Common across items", "Equal (=1)", 
           "Least complex", "Rasch family", "Rating scales"),
    
    PCM = c("Adjacent categories", "Item-specific", "Equal (=1)", 
           "Moderate", "Rasch family", "Achievement tests"),
    
    GPCM = c("Adjacent categories", "Item-specific", "Item-specific", 
            "Complex", "Unidimensional", "Mixed item types")
  )
  
  print(model_comparison)
  
  return(list(
    rsm_data = rsm_responses,
    pcm_data = pcm_responses,
    true_parameters = list(
      theta = theta,
      item_difficulties = item_difficulties,
      common_thresholds = common_thresholds,
      item_thresholds = item_thresholds
    ),
    usage_analysis = list(rsm = rsm_usage, pcm = pcm_usage),
    fit_comparison = if(exists("fit_comparison")) fit_comparison else NULL
  ))
}

# Run RSM and PCM demonstration
rsm_pcm_results <- demonstrate_rsm_pcm_models()
## === RSM VƏ PCM MODELS DEMONSTRASIYASI ===
## Model Configurations:
## RSM - Common thresholds: -1.2, 0, 1.2 
## PCM - Item-specific thresholds (first 3 items):
##   Item 1 : -0.29, 0.58, 0.66 
##   Item 2 : -0.82, -0.59, -0.1 
##   Item 3 : -0.56, -0.11, 1.54 
## 
## Response Data Generated:
## RSM data range: 0 3 
## PCM data range: 0 3 
## 
## === DESCRIPTIVE ANALYSIS ===
## Category Usage - RSM :
## Overall: 1940 / 2535 / 2723 / 2402 ( 20.2% / 26.4% / 28.4% / 25 %)
## Item-level usage (first 4 items):
##    Item Cat0 Cat1 Cat2 Cat3 Mean   SD
## 0     1   25  136  387  652 2.39 0.77
## 01    2  276  410  346  168 1.34 0.98
## 02    3  235  417  350  198 1.43 0.98
## 03    4  103  299  434  364 1.88 0.94
## 
## Category Usage - PCM :
## Overall: 2368 / 2193 / 2138 / 2901 ( 24.7% / 22.8% / 22.3% / 30.2 %)
## Item-level usage (first 4 items):
##    Item Cat0 Cat1 Cat2 Cat3 Mean   SD
## 0     1   50  144  256  750 2.42 0.86
## 01    2  279  260  274  387 1.64 1.16
## 02    3  377  306  365  152 1.24 1.03
## 03    4  175  404  306  315 1.63 1.02
## 
## === MODEL FITTING VƏ COMPARISON ===
## Fitting models with eRm package...
## RSM fitted successfully
## PCM fitted successfully
## 
## Fitting models with mirt package...
## RSM fitted with mirt
## PCM fitted with mirt
## GPCM fitted with mirt
## 
## Model Fit Comparison:
##   Model      AIC      BIC    LogLik Delta_AIC Delta_BIC
## 1   RSM 20980.31 21036.31 -10479.16    0.0000    0.0000
## 2   PCM 21504.92 21632.17 -10727.46  524.6035  595.8645
## 3  GPCM 21514.12 21677.00 -10725.06  533.8038  640.6955
## 
## Parameter Comparison (first 3 items):
## RSM item parameters (sample):
##       a1     b1     b2     b3      c
## Item1  1 -3.008 -1.719 -0.482  0.000
## Item2  1 -3.008 -1.719 -0.482 -2.018
## Item3  1 -3.008 -1.719 -0.482 -1.866
## 
## GPCM item parameters (sample):
##          a1 ak0 ak1 ak2 ak3 d0    d1    d2     d3
## Item1 0.958   0   1   2   3  0 2.201 3.361  4.407
## Item2 0.936   0   1   2   3  0 0.510 0.634  0.545
## Item3 0.992   0   1   2   3  0 0.227 0.242 -1.374
## 
## === MODEL CHARACTERISTICS COMPARISON ===
##                Aspect                   GRM                 RSM
## 1    Parameterization     Cumulative logits Adjacent categories
## 2 Threshold Structure Item-specific ordered Common across items
## 3      Discrimination         Item-specific          Equal (=1)
## 4          Complexity          Most complex       Least complex
## 5         Assumptions        Unidimensional        Rasch family
## 6            Best For         Likert scales       Rating scales
##                   PCM                GPCM
## 1 Adjacent categories Adjacent categories
## 2       Item-specific       Item-specific
## 3          Equal (=1)       Item-specific
## 4            Moderate             Complex
## 5        Rasch family      Unidimensional
## 6   Achievement tests    Mixed item types
cat("\n=== MODEL SELECTION GUIDELINES ===\n")
## 
## === MODEL SELECTION GUIDELINES ===
# Model selection guidelines
selection_guidelines <- c(
  "**MODEL SELECTION CRITERIA:**",
  "",
  "1. **Theoretical Considerations:**",
  "   • GRM: For Likert-type scales with ordered responses",
  "   • RSM: When threshold structure should be common across items",
  "   • PCM: For achievement tests with item-specific difficulty progression",
  "   • GPCM: When items differ in discrimination",
  "",
  "2. **Data Characteristics:**",
  "   • Sample size: RSM/PCM need less data than GRM/GPCM",
  "   • Category usage: All models require adequate category frequencies",
  "   • Item heterogeneity: More varied items favor GRM/GPCM",
  "",
  "3. **Statistical Criteria:**",
  "   • Model fit: Chi-square, AIC, BIC comparisons",
  "   • Parameter stability: Reasonable estimates and standard errors",
  "   • Residual analysis: Person and item fit statistics",
  "",
  "4. **Practical Considerations:**",
  "   • Interpretability: Simpler models often preferred",
  "   • Software availability: Not all packages support all models",
  "   • Reporting requirements: Some contexts prefer specific approaches"
)

for(guideline in selection_guidelines) {
  cat(guideline, "\n")
}
## **MODEL SELECTION CRITERIA:** 
##  
## 1. **Theoretical Considerations:** 
##    • GRM: For Likert-type scales with ordered responses 
##    • RSM: When threshold structure should be common across items 
##    • PCM: For achievement tests with item-specific difficulty progression 
##    • GPCM: When items differ in discrimination 
##  
## 2. **Data Characteristics:** 
##    • Sample size: RSM/PCM need less data than GRM/GPCM 
##    • Category usage: All models require adequate category frequencies 
##    • Item heterogeneity: More varied items favor GRM/GPCM 
##  
## 3. **Statistical Criteria:** 
##    • Model fit: Chi-square, AIC, BIC comparisons 
##    • Parameter stability: Reasonable estimates and standard errors 
##    • Residual analysis: Person and item fit statistics 
##  
## 4. **Practical Considerations:** 
##    • Interpretability: Simpler models often preferred 
##    • Software availability: Not all packages support all models 
##    • Reporting requirements: Some contexts prefer specific approaches
cat("\n=== POLYTOMOUS MODEL RECOMMENDATIONS ===\n")
## 
## === POLYTOMOUS MODEL RECOMMENDATIONS ===
# Recommendations table
model_recommendations <- data.frame(
  Application = c("Likert Scales", "Rating Scales", "Achievement Tests", 
                 "Mixed Item Types", "Attitude Surveys", "Performance Assessments"),
  
  Primary_Choice = c("GRM", "RSM", "PCM/GPCM", "GPCM", "GRM", "GPCM"),
  
  Alternative = c("GPCM", "GRM", "RSM", "GRM", "RSM", "GRM"),
  
  Key_Consideration = c("Response process", "Common structure", "Item difficulty", 
                       "Item variation", "Response style", "Scoring method"),
  
  Sample_Size = c("Large", "Moderate", "Moderate", "Large", "Moderate", "Large"),
  
  Typical_Categories = c("3-7", "3-5", "2-4", "Variable", "3-7", "3-6")
)

print(model_recommendations)
##               Application Primary_Choice Alternative Key_Consideration
## 1           Likert Scales            GRM        GPCM  Response process
## 2           Rating Scales            RSM         GRM  Common structure
## 3       Achievement Tests       PCM/GPCM         RSM   Item difficulty
## 4        Mixed Item Types           GPCM         GRM    Item variation
## 5        Attitude Surveys            GRM         RSM    Response style
## 6 Performance Assessments           GPCM         GRM    Scoring method
##   Sample_Size Typical_Categories
## 1       Large                3-7
## 2    Moderate                3-5
## 3    Moderate                2-4
## 4       Large           Variable
## 5    Moderate                3-7
## 6       Large                3-6

6 Model Estimation və Practical Applications

6.1 GRM Implementation və Real Data Analysis

cat("=== MODEL ESTIMATION VƏ PRACTICAL APPLICATIONS ===\n\n")
## === MODEL ESTIMATION VƏ PRACTICAL APPLICATIONS ===
cat("Practical GRM Implementation: Real data ilə GRM tətbiqi\n")
## Practical GRM Implementation: Real data ilə GRM tətbiqi
cat("• Parameter estimation strategies\n")
## • Parameter estimation strategies
cat("• Model fit assessment\n")
## • Model fit assessment
cat("• Diagnostic procedures\n")
## • Diagnostic procedures
cat("• Score interpretation və reporting\n\n")
## • Score interpretation və reporting
# Comprehensive practical applications demonstration
demonstrate_practical_applications <- function() {
  
  cat("=== PRACTICAL APPLICATIONS DEMONSTRASIYASI ===\n")
  
  # Create realistic survey data
  set.seed(1901)
  
  # Simulate a psychological scale (e.g., life satisfaction)
  n_items <- 12
  n_respondents <- 800
  n_categories <- 7  # 1-7 Likert scale
  
  # Generate respondent abilities (latent trait: life satisfaction)
  life_satisfaction <- rnorm(n_respondents, 0, 1)
  
  # Create items with different characteristics
  item_params <- data.frame(
    item = 1:n_items,
    discrimination = c(1.2, 1.8, 2.1, 1.5, 0.9, 1.7, 2.3, 1.4, 1.9, 1.1, 1.6, 2.0),
    content_area = rep(c("Overall", "Family", "Work", "Health"), each = 3)
  )
  
  # Generate thresholds (transform to 1-7 scale thinking)
  true_thresholds <- matrix(NA, n_items, n_categories - 1)
  
  for(i in 1:n_items) {
    # Generate ordered thresholds, then adjust for different response tendencies
    base_thresholds <- seq(-2.5, 2.5, length.out = n_categories - 1)
    
    # Add some variation
    noise <- rnorm(n_categories - 1, 0, 0.2)
    item_thresholds <- sort(base_thresholds + noise)
    
    true_thresholds[i, ] <- item_thresholds
  }
  
  # Generate responses (1-7 scale)
  responses <- matrix(NA, n_respondents, n_items)
  
  for(i in 1:n_respondents) {
    for(j in 1:n_items) {
      
      # Convert to 0-6 for calculation, then add 1
      cat_probs <- grm_probability(life_satisfaction[i], 
                                  item_params$discrimination[j], 
                                  true_thresholds[j, ])
      
      response_0to6 <- sample(0:(n_categories-1), 1, prob = cat_probs)
      responses[i, j] <- response_0to6 + 1  # Convert to 1-7 scale
    }
  }
  
  colnames(responses) <- paste0("LS", sprintf("%02d", 1:n_items))
  
  cat("Life Satisfaction Scale Simulation:\n")
  cat("Respondents:", n_respondents, "\n")
  cat("Items:", n_items, "\n")
  cat("Scale range:", range(responses), "\n")
  cat("Content areas:", paste(unique(item_params$content_area), collapse = ", "), "\n\n")
  
  # Descriptive analysis
  cat("=== DESCRIPTIVE ANALYSIS ===\n")
  
  # Item descriptives
  item_descriptives <- data.frame(
    Item = colnames(responses),
    Mean = round(apply(responses, 2, mean), 2),
    SD = round(apply(responses, 2, sd), 2),
    Min = apply(responses, 2, min),
    Max = apply(responses, 2, max),
    Skewness = round(apply(responses, 2, function(x) {
      n <- length(x)
      mean_x <- mean(x)
      sd_x <- sd(x)
      skew <- (n / ((n - 1) * (n - 2))) * sum(((x - mean_x) / sd_x)^3)
      return(skew)
    }), 2)
  )
  
  cat("Item Descriptives (first 6 items):\n")
  print(item_descriptives[1:6, ])
  
  # Category frequencies
  cat("\nCategory Usage Analysis:\n")
  overall_freq <- table(factor(as.vector(responses), levels = 1:7))
  overall_prop <- prop.table(overall_freq)
  
  cat("Overall category usage:\n")
  for(i in 1:7) {
    cat("Category", i, ":", overall_freq[i], 
        "(", round(overall_prop[i] * 100, 1), "%)\n")
  }
  
  # Reliability analysis
  cat("\n=== RELIABILITY ANALYSIS ===\n")
  
  if(require(psych, quietly = TRUE)) {
    tryCatch({
      reliability_results <- alpha(responses)
      cat("Cronbach's Alpha:", round(reliability_results$total$raw_alpha, 3), "\n")
      cat("Standardized Alpha:", round(reliability_results$total$std.alpha, 3), "\n")
      
      # Item-total correlations
      item_total_corr <- cor(responses, rowSums(responses))
      cat("Item-total correlations (first 6):", 
          paste(round(item_total_corr[1:6], 3), collapse = ", "), "\n")
      
    }, error = function(e) {
      cat("Reliability analysis failed:", e$message, "\n")
    })
  }
  
  # Fit GRM using mirt
  cat("\n=== GRM MODEL FITTING ===\n")
  
  if(require(mirt, quietly = TRUE)) {
    
    tryCatch({
      
      # Convert to 0-6 scale for mirt (mirt expects 0-based categories)
      responses_0to6 <- responses - 1
      
      # Fit GRM
      cat("Fitting Graded Response Model...\n")
      grm_fit <- mirt(responses_0to6, 1, itemtype = "graded", verbose = FALSE)
      
      cat("GRM fitted successfully\n")
      
      # Extract fit statistics
      fit_stats <- data.frame(
        Statistic = c("Log-likelihood", "AIC", "BIC", "Number of Parameters"),
        Value = c(
          round(extract.mirt(grm_fit, "logLik"), 2),
          round(extract.mirt(grm_fit, "AIC"), 2),
          round(extract.mirt(grm_fit, "BIC"), 2),
          extract.mirt(grm_fit, "nfact") * n_items + sum(apply(responses_0to6, 2, max))
        )
      )
      
      cat("\nModel Fit Statistics:\n")
      print(fit_stats)
      
      # Extract parameter estimates
      grm_coefs <- coef(grm_fit, simplify = TRUE)
      
      # Item parameters
      item_estimates <- data.frame(
        Item = colnames(responses),
        Discrimination = round(grm_coefs$items[, "a1"], 3),
        Threshold_1 = round(grm_coefs$items[, "d1"] / (-grm_coefs$items[, "a1"]), 3),
        Threshold_2 = round(grm_coefs$items[, "d2"] / (-grm_coefs$items[, "a1"]), 3),
        Threshold_3 = round(grm_coefs$items[, "d3"] / (-grm_coefs$items[, "a1"]), 3),
        Content = item_params$content_area
      )
      
      cat("\nItem Parameter Estimates (first 6 items):\n")
      print(item_estimates[1:6, ])
      
      # Parameter recovery analysis
      cat("\n=== PARAMETER RECOVERY ANALYSIS ===\n")
      
      estimated_disc <- grm_coefs$items[, "a1"]
      true_disc <- item_params$discrimination
      
      disc_correlation <- cor(estimated_disc, true_disc)
      disc_rmse <- sqrt(mean((estimated_disc - true_disc)^2))
      
      cat("Discrimination parameter recovery:\n")
      cat("Correlation:", round(disc_correlation, 3), "\n")
      cat("RMSE:", round(disc_rmse, 3), "\n")
      
      # Threshold recovery (first threshold only for simplicity)
      estimated_thresh1 <- grm_coefs$items[, "d1"] / (-grm_coefs$items[, "a1"])
      true_thresh1 <- true_thresholds[, 1]
      
      thresh_correlation <- cor(estimated_thresh1, true_thresh1)
      thresh_rmse <- sqrt(mean((estimated_thresh1 - true_thresh1)^2))
      
      cat("First threshold recovery:\n")
      cat("Correlation:", round(thresh_correlation, 3), "\n")
      cat("RMSE:", round(thresh_rmse, 3), "\n")
      
      # Person parameter estimates (ability scores)
      cat("\n=== PERSON PARAMETER ESTIMATION ===\n")
      
      person_scores <- fscores(grm_fit, method = "EAP")
      
      score_correlation <- cor(person_scores[, 1], life_satisfaction)
      score_rmse <- sqrt(mean((person_scores[, 1] - life_satisfaction)^2))
      
      cat("Person ability estimation:\n")
      cat("Correlation with true θ:", round(score_correlation, 3), "\n")
      cat("RMSE:", round(score_rmse, 3), "\n")
      
      # Score distribution
      cat("Score distribution:\n")
      cat("Mean:", round(mean(person_scores[, 1]), 3), "\n")
      cat("SD:", round(sd(person_scores[, 1]), 3), "\n")
      cat("Range:", paste(round(range(person_scores[, 1]), 2), collapse = " to "), "\n")
      
      # Model diagnostics
      cat("\n=== MODEL DIAGNOSTICS ===\n")
      
      # Item fit
      item_fit <- itemfit(grm_fit, method = "S_X2")
      
      if(!is.null(item_fit)) {
        significant_misfit <- sum(item_fit$p.S_X2 < 0.01, na.rm = TRUE)
        cat("Items with significant misfit (p < 0.01):", significant_misfit, "out of", n_items, "\n")
        
        if(significant_misfit > 0) {
          misfit_items <- which(item_fit$p.S_X2 < 0.01)
          cat("Misfitting items:", paste(colnames(responses)[misfit_items], collapse = ", "), "\n")
        }
      }
      
      # Person fit
      person_fit <- personfit(grm_fit)
      
      if(!is.null(person_fit)) {
        significant_person_misfit <- sum(person_fit$p.Zh < 0.01, na.rm = TRUE)
        cat("Persons with significant misfit (p < 0.01):", significant_person_misfit, 
            "out of", n_respondents, "\n")
      }
      
      # Create practical interpretation
      cat("\n=== PRACTICAL INTERPRETATION ===\n")
      
      # Score conversion table
      theta_levels <- seq(-2, 2, 0.5)
      expected_scores <- numeric(length(theta_levels))
      
      for(i in 1:length(theta_levels)) {
        item_expected_scores <- numeric(n_items)
        
        for(j in 1:n_items) {
          cat_probs <- grm_probability(theta_levels[i], 
                                      estimated_disc[j], 
                                      grm_coefs$items[j, paste0("d", 1:6)] / (-estimated_disc[j]))
          item_expected_scores[j] <- sum(cat_probs * (0:6))
        }
        
        expected_scores[i] <- sum(item_expected_scores)
      }
      
      score_conversion <- data.frame(
        Theta_Level = theta_levels,
        Expected_Total_Score = round(expected_scores, 1),
        Interpretation = c("Very Low", "Low", "Below Average", "Average", 
                          "Above Average", "High", "Very High", "Extremely High", "Maximum")[1:length(theta_levels)]
      )
      
      cat("Score Interpretation Table:\n")
      print(score_conversion)
      
      return(list(
        responses = responses,
        grm_fit = grm_fit,
        item_estimates = item_estimates,
        person_scores = person_scores,
        fit_stats = fit_stats,
        recovery = list(
          disc_corr = disc_correlation,
          thresh_corr = thresh_correlation,
          score_corr = score_correlation
        ),
        score_conversion = score_conversion
      ))
      
    }, error = function(e) {
      cat("GRM fitting failed:", e$message, "\n")
      return(NULL)
    })
    
  } else {
    cat("mirt package not available\n")
    return(NULL)
  }
}

# Run practical applications demonstration
practical_results <- demonstrate_practical_applications()
## === PRACTICAL APPLICATIONS DEMONSTRASIYASI ===
## Life Satisfaction Scale Simulation:
## Respondents: 800 
## Items: 12 
## Scale range: 1 7 
## Content areas: Overall, Family, Work, Health 
## 
## === DESCRIPTIVE ANALYSIS ===
## Item Descriptives (first 6 items):
##      Item Mean   SD Min Max Skewness
## LS01 LS01 3.90 1.58   1   7     0.07
## LS02 LS02 3.96 1.45   1   7    -0.07
## LS03 LS03 3.98 1.27   1   7     0.22
## LS04 LS04 3.98 1.42   1   7     0.09
## LS05 LS05 4.00 1.84   1   7     0.00
## LS06 LS06 3.82 1.34   1   7     0.06
## 
## Category Usage Analysis:
## Overall category usage:
## Category 1 : 601 ( 6.3 %)
## Category 2 : 860 ( 9 %)
## Category 3 : 1998 ( 20.8 %)
## Category 4 : 2781 ( 29 %)
## Category 5 : 1994 ( 20.8 %)
## Category 6 : 865 ( 9 %)
## Category 7 : 501 ( 5.2 %)
## 
## === RELIABILITY ANALYSIS ===
## Cronbach's Alpha: 0.883 
## Standardized Alpha: 0.89 
## Item-total correlations (first 6): 0.594, 0.715, 0.743, 0.667, 0.555, 0.669 
## 
## === GRM MODEL FITTING ===
## Fitting Graded Response Model...
## GRM fitted successfully
## 
## Model Fit Statistics:
##              Statistic     Value
## 1       Log-likelihood -14688.47
## 2                  AIC  29544.94
## 3                  BIC  29938.45
## 4 Number of Parameters     84.00
## 
## Item Parameter Estimates (first 6 items):
##      Item Discrimination Threshold_1 Threshold_2 Threshold_3 Content
## LS01 LS01          1.155      -2.346      -1.740      -0.600 Overall
## LS02 LS02          1.757      -2.086      -1.729      -0.368 Overall
## LS03 LS03          2.151      -2.411      -1.802      -0.451 Overall
## LS04 LS04          1.569      -2.581      -1.636      -0.415  Family
## LS05 LS05          0.927      -2.516      -1.228      -0.878  Family
## LS06 LS06          1.612      -2.381      -1.708      -0.289  Family
## 
## === PARAMETER RECOVERY ANALYSIS ===
## Discrimination parameter recovery:
## Correlation: 0.976 
## RMSE: 0.091 
## First threshold recovery:
## Correlation: 0.754 
## RMSE: 0.167 
## 
## === PERSON PARAMETER ESTIMATION ===
## Person ability estimation:
## Correlation with true θ: 0.952 
## RMSE: 0.305 
## Score distribution:
## Mean: 0 
## SD: 0.952 
## Range: -2.82 to 2.51 
## 
## === MODEL DIAGNOSTICS ===
## GRM fitting failed: unused argument (method = "S_X2")
cat("\n=== PRACTICAL IMPLEMENTATION CHECKLIST ===\n")
## 
## === PRACTICAL IMPLEMENTATION CHECKLIST ===
# Implementation checklist
implementation_checklist <- c(
  "**PRE-ANALYSIS CHECKS:**",
  "☐ Adequate sample size (200+ respondents minimum)",
  "☐ Sufficient category usage (5+ responses per category)",
  "☐ Data quality (missing data patterns, outliers)",
  "☐ Assumption verification (unidimensionality, local independence)",
  "",
  "**MODEL FITTING:**",
  "☐ Software selection and parameterization",
  "☐ Starting value specification",
  "☐ Convergence monitoring",
  "☐ Parameter estimate reasonableness",
  "",
  "**MODEL EVALUATION:**",
  "☐ Overall model fit assessment",
  "☐ Item-level fit diagnostics",
  "☐ Person-level fit analysis",
  "☐ Parameter stability across subgroups",
  "",
  "**INTERPRETATION AND REPORTING:**",
  "☐ Score conversion and interpretation guidelines",
  "☐ Reliability and precision estimates",
  "☐ Limitation and assumption documentation",
  "☐ Practical significance assessment"
)

for(check in implementation_checklist) {
  cat(check, "\n")
}
## **PRE-ANALYSIS CHECKS:** 
## ☐ Adequate sample size (200+ respondents minimum) 
## ☐ Sufficient category usage (5+ responses per category) 
## ☐ Data quality (missing data patterns, outliers) 
## ☐ Assumption verification (unidimensionality, local independence) 
##  
## **MODEL FITTING:** 
## ☐ Software selection and parameterization 
## ☐ Starting value specification 
## ☐ Convergence monitoring 
## ☐ Parameter estimate reasonableness 
##  
## **MODEL EVALUATION:** 
## ☐ Overall model fit assessment 
## ☐ Item-level fit diagnostics 
## ☐ Person-level fit analysis 
## ☐ Parameter stability across subgroups 
##  
## **INTERPRETATION AND REPORTING:** 
## ☐ Score conversion and interpretation guidelines 
## ☐ Reliability and precision estimates 
## ☐ Limitation and assumption documentation 
## ☐ Practical significance assessment
cat("\n=== COURSE COMPLETION ===\n")
## 
## === COURSE COMPLETION ===
cat("Bu dərsi tamamlayaraq, siz Graded Response Model və polytomous IRT-nin əsas prinsiplərini öyrəndiniz.\n")
## Bu dərsi tamamlayaraq, siz Graded Response Model və polytomous IRT-nin əsas prinsiplərini öyrəndiniz.
cat("İndi Likert scales və rating data-nı sophisticated IRT methods ilə analiz edə bilərsiniz.\n\n")
## İndi Likert scales və rating data-nı sophisticated IRT methods ilə analiz edə bilərsiniz.
cat("Next Steps for Polytomous IRT Mastery:\n")
## Next Steps for Polytomous IRT Mastery:
cat("• Practice with ltm, mirt, and eRm packages\n")
## • Practice with ltm, mirt, and eRm packages
cat("• Explore nominal response models for unordered categories\n")
## • Explore nominal response models for unordered categories
cat("• Study multidimensional polytomous IRT models\n")
## • Study multidimensional polytomous IRT models
cat("• Learn about mixture polytomous IRT models\n")
## • Learn about mixture polytomous IRT models
cat("• Investigate computerized adaptive testing with polytomous items\n")
## • Investigate computerized adaptive testing with polytomous items
cat("• Apply polytomous IRT to real survey and assessment data\n")
## • Apply polytomous IRT to real survey and assessment data
cat("• Compare polytomous IRT with classical methods\n\n")
## • Compare polytomous IRT with classical methods
cat("Polytomous IRT opens new possibilities for sophisticated measurement!\n")
## Polytomous IRT opens new possibilities for sophisticated measurement!
cat("Müvəffəqiyyətlər polytomous IRT expertise journey-nizdə!\n")
## Müvəffəqiyyətlər polytomous IRT expertise journey-nizdə!